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ABSTRACT 

For  advection  schemes  based  on  fluctuation  splitting,  a  design  criterion  of  optimising  the 
time  step  leads  to  linear  schemes  that  coincide  with  those  designed  for  least  truncation  error. 
A  further  stage  of  optimising  the  time  step  using  a  non-linear  positivity  criterion,  leads  to 
considerable  further  gains  in  resolution. 
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1.  Introduction 


This  note  is  a  contribution  to  the  development  of  high  resolution  advection  schemes  for 
use  on  (possibly)  unstructured  tricingular  meshes.  The  seeirch  is  for  a  scheme  that  is  positive 
(admitting  no  new  extrema)  and  as  accurate  as  possible.  Such  schemes  have  an  inherently 
nonUnear  dependence  on  the  data,  and  revert  near  discontinuities  to  an  essentially  first-order 
behavior.  It  is  of  interest  to  know  what  is  the  optimum  positive  first-order  scheme,  so  that 
performance  near  discontinuities  is  not  excessively  degraded. 

The  schemes  considered  in  this  paper  belong  to  a  class  introduced  in  [1],  and  explored 
more  thoroughly  in  [2].  These  papers,  and  also  [3],  should  be  consulted  for  more  background. 

2.  Previous  Work 

Consider  a  trianguleir  grid  of  which  Figure  1  shows  a  generic  element.  On  such  a  grid, 
we  wish  to  solve  the  two-dimensional  advection  equation 


ut  +  a  ■  Vu  =  0, 


where  a  is  a  constant  vector  (a,  6)^. 

For  each  trianguleir  element  T  define  the  “fluctuation” 

(f)^  J  J  Utdxdy. 

By  Equation  (1)  and  Gauss’  theorem 

(f>^  —  ^  u  a  -  dn. 

If  u  is  assumed  to  vary  linearly  over  the  triangle 

f 

where 

.  1 

Ki  =  -a-  Hi, 

and  rii  is  the  inward  (scaled)  normal  to  the  side  opposite  Vi.  Note  that 

=  0 


(1) 


(2) 

(3) 


so  that  either  one  or  two  ki  are  positive.  Geometrically,  this  corresponds  to  the  vector  a 
crossing  either  one  or  two  sides  in  the  inward  direction. 

In  [1 ,2]  it  is  shown  that  part  of  a  positive  update  process  is  to  add,  in  the  case  of  one  inflow 
side  (which  we  call  a  type  I  triangle.  Figure  2a)  a  quantity  At<f>T  to  the  single  downstream 
vertex.  If  this  is  Vi,  then 

=  Siu:  Atf  -h  (■)  (4) 

where  (•)  signifies  the  contributions,  if  any,  coming  from  other  triangles  neighboring  Vi,  and 
Si  is  one-third  the  area  of  those  triangles  having  as  a  vertex. 
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In  the  difficult  case  where  there  are  two  inflow  sides,  say  those  opposite  vertices  i  and  j 
(a  type  II  triangle)  the  proposed  scheme  is 


=  5iuJ‘  +  +  (•)  \ 

=  Sju^  +  +  (-)  I 


(5) 


where,  from  a  conservation  argument. 


ai  +  Oj  =  1. 


It  was  proved  in  [1]  (see  also  [2])  that  no  choice  of  o-j,  aj  depending  only  on  the 
of  the  triangle  and  its  relation  to  the  flow  vector  can  give  a  positive  scheme.  The 
choice,  which  we  call  Si, 

ai  =  ki{ui  -  Uk)l(f>^  I 


geometry 

nonlinear 


(6) 


cij  =  kj{uj  -  Uk)l<i>^  J 

was  advocated  in  [1]  on  the  grounds  that  it  permitted  a  larger  timestep  than  any  other  choice 
(within  a  certain  class).  This  scheme  was  found  to  be  superior  to  regular  finite- volume 
upwinding  in  [2,  3]  and  we  will  point  out  in  the  next  section  that  it  bears  an  interesting 
relationship  to  a  recently  published  scheme  of  Sidilkover  [4]. 


3.  A  Property  of  the  Scheme  Si 


Substituting  (6)  into  (5)  gives  the  simple  update  method 

=  SiU^  -  AiAt(tir  -  njj)  -f  (•)  I 

-  kjAt{u^  -  u'i^)  +  (•)  J 


(7) 


Let  us  see  what  (7)  reduces  to  on  a  square  mesh  of  side  h.  We  assume  that  a,  b  are  both 
positive  with  a  >  b,  and  that  the  squares  are  divided  into  triangles  along  the  diagonal  with 
positive  slope,  as  in  Figure  3a.  Freedom  to  cut  the  squares  along  either  diagonal  greatly 
enhances  the  performance  of  the  schemes  (R.  Struijs,  H.  Deconinck,  private  communication, 
von  Karman  Institute). 

Each  point  then  receives  contributions  from  three  triangles.  A,  B,  C,  in  which  the  fluctu¬ 
ations  are  ^ 

<f)^  =  --[(a  -  b){uo  -  U3)h  -f  b{ui  -  U3)h] 

=  -^[a(uo  -  U3)h  -I-  b{u3  -  Ui)h] 

<i)^  =  -^[(a  -  6)(u5  -  u^)h  ->r  bjup  -  Ui)h]. 

In  each  case,  the  fluctuations  have  been  split  as  they  would  be  by  SI,  and  those  terms 
contributing  to  the  update  of  uq  have  been  underlined.  Summing  those  terms,  we  find 

=  ■Uq  -  ~  ^3)  +  ^“3  -  «4)]-  (8) 


2 


This  is  the  scheme  recently  shown  by  Sidilkover  [4]  to  have  the  least  truncation  error 
among  the  linear  monotone  schemes  for  solving  (1).  It  comes  in  eight  different  flavors 
depending  on  the  flow  direction  (Figure  4).  (It  may  be  remarked  that  whether  the  scheme 
is  linear  or  nonlinear  depends  on  the  viewpoint.  Clearly  (7)  has  constant  coefficients,  but 
they  change  as  the  vector  o  moves  into  different  octants.  Written  in  the  form  (5),  (6)  the 
expressions  that  accomplish  this  switching  give  the  appearance  of  nonlinearity.)  Sidilkover 
[4]  calls  (7)  the  JV-scheme  because  it  has  a  narrow  stencil  and,  not  by  coincidence,  yields 
narrow  discontinuities. 

Now  consider  the  squares  divided  along  the  other  family  of  diagonals  as  in  Figure  3b. 
The  fluctuations  in  the  triangles  adjoining  0  are 

(f)^  —  -a{uQ  -  1x3)  -  6(u2  -  U3) 

4)^  =  -a(uo  -  tta)  -  &(uo  -  u^) 

=  -a{ue  -  Us)  -  b{uQ  -  1x5) 

and  scheme  SI  will  credit  the  underlined  terms  to  0.  The  update  formula  is 

-  ~[a{uo  -  U3)  +  b(uo  -  Us)]  (9) 

which  is  just  regular  upwinding.  For  steady  flows,  the  dissipation  coefficient  of  (9)  is 

^|sin5cos^|{|sin5|  +  jcos^l} 


whereas  the  dissipation  coefficient  of  (8)  is 


—  I  sin  0  cos  I  cos  0  —  sin  0\ 


which  is,  on  average,  about  four  times  smaller.  There  seems  no  reason  why  one  should  not 
taJce  advantage  of  this  by  choosing  the  triangulation  appropriate  to  the  problem.  In  the 
non-linear  case,  of  course,  the  choice  might  change  with  time. 


4.  Type  II  Triangles  Revisited 

The  observation  is  made  in  [1,  2]  that  within  each  triangle  T  the  advection  velc  .ty  a  can 
be  replaced  by 

=  a  +  AoJ  (10) 

where  oj  is  orthogonal  to  Vu^  and  A  is  an  arbitrary  parameter.  This  is  tar^^amount  to  saying 
that  only  the  component  of  a  normal  to  the  level  lines  of  the  solution  i'.  relevant  (Figure  5). 
In  any  analysis  relating  to  events  within  the  triangle  T,  a  can  be  rr, 'laced  by  o^.  There  is 
an  intuitive  attraction  to  making  as  small  as  possible,  that  is  taking  it  in  the  direction 
Vu^.  Then  we  have 

T  (  a- j. 

~  \Vu'^  ■  -  vu-  .  Vu^  ■ 
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A  practical  disadvantage  to  (9)  is  that  in  any  region  where  the  solution  is  almost  uniform, 
the  numerator  ajid  denominator  are  both  very  small,  and  numerical  expressions  based  on 
(9)  become  ill-formed.  This  leads  to  slow,  or  halted,  convergence.  In  [2],  the  expedient  was 
taken  of  abandoning  (9)  gradually  in  favor  of  a  in  regions  where  |V„  j  is  small.  Below,  we 
will  argue  that  (9)  is  anyway  not  optimal;  the  improved  version  leads  also  to  better  formed 
expressions. 

We  commence  the  analysis  by  noting  that  if  (10)  is  substituted  into  (3)  we  have 

ki  =  •  Tii 


=  a  •  ni  +  -H.- 

In  the  second  inner  product,  is  orthogonal  to  Vu  and  is  orthogonal  to  where 
is  the  vector  along  the  side  opposite  Vi).  Therefore 

ki  —  a  -  n,  -f-  An  •  (12) 


so  that 

ki  =  k*  +  X(u3  -  U2) 
k^  =  k^  +  A(ui  —  tia)  • 
k3  =  kl  +  X{u2  -  ui)  ^ 


(13) 


Here,  A:J,  /sj,  k^  are  the  original  coefficients  computed  from  (3),  and  equations  (13)  define 
a  one-parameter  family  of  equally  valid  alternatives. 

We  pose  the  problem  of  choosing  the  A:’s  so  as  to  maximize  the  allowable  time-step  of 
the  scheme.  This  is  clearly  a  desirable  objective  in  itself  but  heuristically  it  seems  to  lead  to 
gains  in  accuracy  also.  (In  one  dimension,  this  criterion  would  lead  us  to  classical  upwinding 
and  in  Section  3,  we  saw  that  it  led  to  the  optimum  linear  schemes  in  two  dimensions.)  From 
(7),  the  scheme  is  positive  for  Type  II  triangles,  if 


At  <  —  min 
n 


where  i,j  are  the  downwind  vertices,  and  where  n  is  the  majcimum  number  of  triangles 
meeting  at  a  point.  Therefore,  the  relevant  criterion  is  that 


K  =  max 


(14) 


should  be  positive  but  as  small  as  possible. 

Consider  a  diagram  of  which  the  axes  are  (Aj./S’,),  (Aj/^j).  As  A  varies,  a  straight  line 
locus  will  be  swept  out  whose  slope  is 

A  _  ^«(^i  -  ^fc) 
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Figure  6  shows  such  loci  for  given  k*,k*j  (i.e.,  given  geometry)  for  a  range  of  A  (i.e.,  a 
range  of  data).  For  each  locus,  a  star  marks  the  point  that  minimizes  K. 

If  the  star  falls  on  the  vertical  axis,  the  decision  is  to  update  only  the  point  j\  on  the 
horizontal  axis  we  update  only  the  point  i.  On  the  45°  line  we  have 

k 

\  ^3 

so  that,  from  (7) 

8u^  Ui  —  Uk 

6uJ  Uj  —  Uk 

(where  SuJ  means  the  change  at  vertex^  due  to  this  triangle).  This  corresponds  to  keeping 
the  direction  of  Vu  constant  during  the  update  process.  Figure  6  demonstrates  that  this 
is  the  best  policy  whenever  A  is  negative,  i.e.,  [ui  —  Uk),{uj  —  Uk)  have  the  same  sign.  An 
equivalent  criterion  is  that  the  level  line  of  the  solution  through  Vk  does  not  pass  through 
the  triangle.  If  the  level  line  does  pass  through  the  triangle,  v/e  will  only  update  one  vertex. 
Which  vertex  is  chosen  will  depend  whether  the  locus  in  Figure  6  passes  above  or  below  the 
origin,  and  it  is  easy  to  demonstrate  that  this  depends  on  the  sign  of 
The  update  algorithm  is  in  fact,  in  pseudo- Fortran 

If  (only  k^.gt.O) 
update  Vi 
goto  end 
endif 

If  (only  kk.lt.Q)  then 

If  ((uj  —  Uk)(uj  —  Uk).gt.O)  then 
update  Ui,Uj 

else  if  {{ui  —  Uk)<l>^ .It-O)  then 
update  Ui 
else 

update  Uj 
endif 
endif . 


There  are  actually  twelve  different  paths  through  this  logic,  although  only  six  different 
outcomes.  It  can  be  verified  that  the  updates  at  the  vertices  depend  C°  continuously  on  all 
aspects  of  the  data.  The  only  part  of  the  algorithm  that  is  badly  formed  arises  if  both  Ui,Uj 
are  to  be  updated,  in  accordance  with  (15).  The  explicit  formulae  are  then 


Si{Ui  -  Uk)  +  Sj{uj  -  Uk)^ 


(16a) 


_ Uj  -  Uk _ 

St{Ui  -  Uk)  +  Sj{Uj  -  Uk) 


(16i) 
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and  it  is  possible  for  the  denominator  to  be  small.  However,  we  only  reach  this  situation  if 
{ui  —  Uk){uj  —  Uk)  are  of  the  same  sign,  so  they  must  be  small  individually,  and  <f>T  must  also 
be  small.  Thus,  that  triangle  is  in  equilibrium,  and  can  be  ignored.  In  practice,  we  may  test 
to  see  if  4't  is  less  than  some  moderate  multiple  of  machine  zero. 

We  call  this  the  NN  scheme,  signifying  non-linear  and  narrow. 

5.  Results  and  Discussion 

The  schemes  will  be  demonstrated  on  two  test  problems  shown  in  Figure  7. 

The  first  problem  is  to  solve  the  equation 

Ut+U^  +  ^Uy  =  0 

on  the  unit  square  with  u  =  1  specified  on  the  left  edge  and  u  =  0  on  the  bottom  edge.  This 
is  the  direction  for  which  the  N,  and  NN,  schemes  will  show  least  advantage.  The  second 
problem  is  to  solve  the  equation 


Ut  +  X  Uy  —  yUx  =  0 

on  the  domain  ne[— 1,1],  t/e[0,  Ij.  The  characteristic  lines  here  are  circles  centered  on  the 
origin. 

On  all  inflow  boundaries  we  specify  u  =  0,  except  that  on  y  =  0,  we  set  u  =  1  for 
ie[— 5,  —'ll.  The  steady  solution  of  this  problem  is  therefore 

u  =  1  if  1  <  9(x*  4-  y*)  <  4 
=  0  otherwise. 

In  particular,  when  we  plot  u{x,  0)  we  expect,  or  at  least  hope,  to  see  the  input  data  for 
X  <  0  mirrored  in  the  solution  for  x  >  0.  An  interesting  feature  of  examining  the  solution  in 
that  way  is  that  the  dissipation  error  has  been  effectively  integrated  over  all  flow  directions. 

Figure  8(a)  contains  results  from  the  regular  upwind  scheme  for  Problem  1  on  a  mesh  of 
31  X  31  nodes.  The  solution  is  very  diffused,  but  would  be  even  worse  at  a  direction  of  45°. 
Figure  8(b)  is  obtained  from  the  SI  scheme,  but  the  mesh  is  triangulated  in  no  preferential 
direction  (Figure  7(c)).  The  results  are  very  slightly  improved,  although  this  cannot  be 
determined  from  the  plots. 

A  very  noticeable  improvement  comes  from  deploying  the  Si  scheme  on  a  mesh  trian¬ 
gulated  along  the  upgoing  diagonals.  The  results  are  shown  in  Figure  8(c),  and  precisely 
the  same  results  are  obtained  using  the  N  scheme  on  a  square  mesh.  The  discontinuity  is 
resolved  over  about  half  the  number  of  mesh  points.  It  is  worth  remarking  that  this  reduces 
by  a  factor  of  about  16  the  computing  cost  of  obtaining  a  solution  with  given  resolution, 
because  half  as  many  points  are  needed  in  each  direction,  half  as  much  time  is  needed  to 
expel  transients,  and  twice  the  timestep  can  be  taken. 

In  Figure  8(d)  we  see  results  from  the  nonlinear  scheme  NN.  Even  on  the  isotropic  mesh, 
these  are  the  best  so  far,  but  by  choosing  the  correct  triangulation,  the  further  improvement 
of  Figure  8(e)  is  reached.  Since  this  improves  the  resolution  by  a  factor  of  4,  the  cost  at 
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given  resolution  is  reduced  by  about  a  factor  of  64.  Recall  that  the  results  from  the  optimum 
scheme  are  actually  worse  for  this  flow  angle.  It  is  gratifying  that  the  convergence  history 
is  scarcely  affected  by  any  of  these  improvements.  The  residual  stays  almost  constant  for 
about  80  timesteps,  by  which  time  the  initial  guess  -u  =  0  has  been  purged  by  the  wave 
crossing  the  domain;  after  a  further  70  to  80  timesteps,  the  solution  is  effectively  converged. 
All  tests  were  made  at  a  Courant  number  of  0.5,  but  in  fact  the  optimum  scheme  could  have 
been  run  faster. 

In  Figure  9(a),  the  results  are  shown  of  the  rotating  advection  problem  on  a  mesh  of  61 
cells  by  31  cells.  Only  the  results  for  an  optimized  mesh  are  shown;  the  squares  are  cut  by 
their  upgoing  diagonals  for  i  <  0,  by  their  downgoing  diagonals  for  i  >  0.  The  square  pulse 
is  preserved  quite  remarkably  for  a  first-order  scheme.  Indeed,  one  may  suspect  that  the  step 
function  is  some  kind  of  eigenmode  for  the  scheme,  and  that  all  data  will  be  squared-off,  as 
happens  for  certain  one- dimensional  limiter  schemes.  To  check  this  point,  a  narrow  Gaussian 
pulse  was  substituted  as  data  and  was  again  preserved  rather  well  (Figure  9(b)). 

To  make  some  kind  of  comparison  with  regular  upwinding,  the  calculation  in  Figure 
9(c)  was  performed.  To  give  the  regular  scheme  a  chance,  four  times  as  many  mesh  points 
were  used  in  each  direction,  i.e.,  241  x  121.  Even  so,  a  very  inferior  result  is  achieved.  In 
comparison  with  the  straight  advection  test,  the  optimum  schemes  show  in  a  better  light  and 
the  regular  scheme  in  a  worse  light.  Their  resolving  powers  differ  here  by  a  factor  of  about 
10,  and  the  cost  of  achieving  given  resolution  changes  by  a  factor  of  about  one  thousa.nd. 

6.  Concluding  Remarks 

The  schemes  derived  here  have  features  with  no  close  one- dimensional  counterparts.  They 
actively  select  the  stencil  to  be  used,  and  the  connectivity  of  the  mesh,  and  depend  non- 
linearly  on  the  data,  all  in  the  search  for  an  optimal  first-order  scheme.  They  are  of  course 
to  be  viewed  as  the  first  phase  in  a  quest  for  higher-order  schemes  that  may  exploit  the  same 
features.  Nevertheless,  their  practical  use  in  simulations  of  mixing  and  turbulent  flows  may 
not  be  out  of  the  question.  The  resolution  that  they  achieve  compares  favorably  with  that 
of  formally  higher-order  schemes,  and  because  all  their  operations  are  highly  localized,  the 
schemes  are  very  suitable  for  parallel  implementation. 
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Figure  3.  Square  cells  divided  either  correctly  (a)  or  incorrectly  (b)  for  flow  with 
positive. 


,  b  both 


Figure  4.  In  the  N-scheme,  point  0  is  updated  using  whichever  trieuigle  contains  the  incoming 
characteristic. 


(b)  Problem  2 


(c)  Section  of  'isotropic'  mesh 


Figure  7.  The  Test  Problems. 
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Figure  8(a)  Problem  1  -  regular  upwinding 
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